%% Data Postprocessing
% Deals with the processing of the run data in various ways

function ProcessData(configs,gui,data)
% Load data and stuff
global globals
if strcmp(configs.MultiFileAnalysis,'None')
if ishandle(12);set(12,'Units','centimeter'); cornerpos = get(12,'Position'); close(12); end;
cfg     = data{1}.cfg;


if isfield(cfg, 'Dataset')
    load(['Data\' cfg.Dataset]);
else
    load('Data\exvars 1971-2100 Original Dataset.mat');
end

%% Pick the best fits
if configs.AllRuns
    configs.Runs = 1:length(data);
end

if configs.Selectbestrun 
    errornames = fieldnames(data{1}.depvars.errors);
    selection  = SelectionMenu(errornames,'Select error indicators to compare');
    errorindic = errornames(selection);
    errorindic = errorindic{1};
    cfg.bestruns = zeros([1, length(cfg.allregions)]);
    errors = zeros([length(cfg.allregions) length(data)]);
        for n = 1 : length(data)
            errors(:,n) = data{n}.depvars.errors.(errorindic);
        end
            [cfg.minimum cfg.bestruns] = min(errors,[],2);

datalength = length(data);
vartypes   = {'depvars', 'indepvars'};
for n = 1:2
   eval(['config.fields = fieldnames( data{1}.' vartypes{n} ');']);
   if n == 1; config.fields = config.fields(~strcmp(config.fields,'errors')); config.fields(end+1:end+length(fieldnames(data{1}.depvars.errors)),1) = cellfun(@(x) ['errors.' x],fieldnames(data{1}.depvars.errors), 'UniformOutput',false); end
   for m = 1:length(config.fields)
       eval(['[matrixsize(1) matrixsize(2) matrixsize(3)] = size(data{1}.' vartypes{n} '.' config.fields{m} ');']);
       for y = 1:size(cfg.bestruns,2)
       eval(['data{datalength+y}.' vartypes{n} '.' config.fields{m} ' = zeros([matrixsize]);']);
       if eval(['size(data{1}.' vartypes{n} '.' config.fields{m} ',1)>1;']);
           for r = 1:length(cfg.allregions)
                eval(['data{datalength+y}.' vartypes{n} '.' config.fields{m} '(r,:,:)=data{cfg.bestruns(r,y)}.' vartypes{n} '.' config.fields{m} '(r,:,:);']);
           end

       end
       end
       disp(config.fields{m});
   end
end

indepvars.UseModes              = exvars.TravelDemand(:,35,:)~=0; % Dump HS trains for regions that don't use it
if cfg.NoUnMotorized || strcmp(cfg.Model,'ModelPoles');    if length(cfg.modes)==4;        indepvars.UseModes(:,:,strcmp(cfg.modes,'NM')) = false;    else;       indepvars.UseModes(:,:,strcmp(cfg.modes,'Walk')|strcmp(cfg.modes,'Bike')) = false;    end;end
indepvars.TTB                   = ones([length(cfg.Regions),length(cfg.Time)])*1.2*365; % hr/yr, seems to fit the data better than increasing with time
indepvars.MinTimeWeight         = ones([length(cfg.Regions),1])*0.1*cfg.RestrictTWRange; 
data{end}.indepvars.UseModes = logical(indepvars.UseModes);
configs.Runs                 = length(data)-size(cfg.bestruns,2)+1:length(data);
data{end}.cfg                = cfg;
end
if ~strcmp(configs.SingleRunMultiRegion,'None')% Used to skip this analysis tools
    r = configs.Regions(1);
    for run = configs.Runs
        %% Transform the regions if desired
        if run==configs.Runs(1)
        originalexvars = exvars;
        end
        [data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,AllRegions] = TransformModesandRegions(data{run}.cfg,originalexvars,data{run}.depvars,data{run}.indepvars,configs.RegionSet);
        [data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars]            = TransformModesandRegions(data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,configs.ModeSet);   
        if configs.AllRegions == 1; configs.Regions = AllRegions; end            
                
                if ishandle(12);set(12,'Units','centimeter'); cornerpos = get(12,'Position'); close(4);end;
                globals.fid = figure(12);set(12,'Units','centimeter')
                if ~exist('cornerpos','var'); cornerpos = get(globals.fid,'Position'); end
        
            switch configs.SingleRunMultiRegion
            
            case {'Barplot indepvars'}
                fields             = fieldnames(data{run}.indepvars);
                selection          = SelectionMenu(fields,'Dimensions to plot');
                selectedfields     = fields(selection);
                
                for n = 1:sum(selection)
                   
                    if size(data{run}.indepvars.(selectedfields{n}),3) >1;
                        modetext = [', mode ' cfg.modes{m}]; mmmm= find(data{run}.indepvars.UseModes(1,:)); else
                        modetext = ''; mmmm = 1; 
                    end
                    for m = mmmm
                    cfg.shortregions = {'NAM','PAO','WEU','EEU','FSU','Afr','CPA','LAM','MEA','PAS','SA','World'};
                    barh(squeeze(data{run}.indepvars.(selectedfields{n})(:,1,m)));
                    set(gca,'yticklabel',cfg.shortregions)
                    xlabel([selectedfields{n}]); ylabel('Region'); title(['Comparison for ' selectedfields{n} modetext]);
                    
                    pubgraph(globals.fid,12,1.5,'white');
                    configs.savename = ['Barh ' selectedfields{n} modetext ' ' datestr(now,'yyyy-mm-dd HH-MM-SS' )];
                    savegraphs(configs);
                    end
                end
                                
            case {'Barplot depvars'}
                data{run}.depvars.TMB        = (sum(data{run}.depvars.TravelDemand.*exvars.TravelPrice./repmat(exvars.Income,[1 1 4]),3));
                data{run}.depvars.TTB        = (sum(data{run}.depvars.TravelDemand./exvars.Speed,3));
                fields             = fieldnames(data{run}.depvars);
                selection          = SelectionMenu(fields,'Dimensions to plot');
                selectedfields     = fields(selection);
                data{run}.depvars.TTB(data{run}.depvars.TTB>1000) = 1000;
                data{run}.depvars.TMB(data{run}.depvars.TMB>0.5) = 0.5;

                
                for n = 1:sum(selection)
                    if size(data{run}.depvars.(selectedfields{n}),3) >1;
                        mmmm= find(data{run}.indepvars.UseModes(1,:)); else
                        mmmm = 1; 
                    end
                    for m = mmmm
                   
                    if size(data{run}.depvars.(selectedfields{n}),3) >1;
                        modetext = [', mode ' cfg.modes{m}]; else
                        modetext = ''; 
                    end
                    
                    cfg.shortregions = {'NAM','PAO','WEU','EEU','FSU','Afr','CPA','LAM','MEA','PAS','SA','World'}; 
                    boxplot(squeeze(data{run}.depvars.(selectedfields{n})(configs.Regions,:,m))','orientation','horizontal','boxstyle','outline','colors','k','outliersize',2,'symbol','kx');
                    h = findobj(gca,'Tag','Box');
                     for j=1:length(h)
                        patch(get(h(j),'XData'),get(h(j),'YData'),'g','FaceAlpha',.5);
                     end
                    set(gca,'yticklabel',cfg.shortregions(configs.Regions),'ytick',[1:12])
                    xlabel([selectedfields{n}]); ylabel('Region'); title(['' selectedfields{n} modetext ]);
                    
                    pubgraph(globals.fid,12,1.5,'white');
                    configs.savename = ['Barh ' cfg.Model selectedfields{n} modetext ' ' datestr(now,'yyyy-mm-dd HH-MM-SS' )];
                    savegraphs(configs);
                    end
                end
                
                case {'Barplot errors'}
                fields             = fieldnames(data{run}.depvars.errors);
                selection          = SelectionMenu(fields,'Dimensions to plot');
                selectedfields     = fields(selection);
                
                for n = 1:sum(selection)
                    for m = 1:size(data{run}.depvars.errors.(selectedfields{n}),3)
                    if size(data{run}.depvars.errors.(selectedfields{n}),3) >1;
                        modetext = [', mode ' cfg.modes{m}]; else
                        modetext = ''; end
                    cfg.shortregions = {'NAM','PAO','WEU','EEU','FSU','Afr','CPA','LAM','MEA','PAS','SA','Wor'};
                    globals.fid = barh(data{run}.depvars.errors.(selectedfields{n})(:,:,m));
                    
                    set(gca,'yticklabel',cfg.shortregions)
                    xlabel([selectedfields{n}]); ylabel('Region'); title(['Comparison for ' selectedfields{n} modetext]);
                    
                    pubgraph(globals.fid,12,1.5,'white');
                    configs.savename = ['Barh ' selectedfields{n} modetext ' ' datestr(now,'yyyy-mm-dd HH-MM-SS' )];
                    savegraphs(configs);
                    end
                end
            end
    end
end


if ~strcmp(configs.SingleRunAnalysis,'None')% Used to skip this analysis tools
    for run = configs.Runs
        %% Transform the regions if desired
        if run==configs.Runs(1)
        originalexvars = exvars;
        end
        
        [data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,AllRegions] = TransformModesandRegions(data{run}.cfg,originalexvars,data{run}.depvars,data{run}.indepvars,configs.RegionSet);
        [data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars]            = TransformModesandRegions(data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,configs.ModeSet);   
        if configs.AllRegions == 1; configs.Regions = AllRegions; end
        
        for r = configs.Regions
        %% Single run analysis tools
            if ishandle(12);set(12,'Units','centimeter'); cornerpos = get(12,'Position'); close(12); end;
            configs.savename = [configs.SingleRunAnalysis ' ' cfg.allregions{r} ' run ' num2str(run) ' ' configs.OutputFile];
            switch configs.SingleRunAnalysis
                case {'Single Plot'}
                    configs.savename = GeneralPlotting(configs.SingleRunAnalysis,'Dont Save',data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,r,1);

                case {'TTB and TMB'}
                    configs.savename = GeneralPlotting(configs.SingleRunAnalysis,'Dont Save',data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,r,1);
                    
                case 'ProjectionArea'
                    % Creates an Area plot, as used int he GCAM publication
                    configs.savename = GeneralPlotting(configs.SingleRunAnalysis,'Dont Save',data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,r,1);

                case {'HistoricPlot','ProjectionPlot','AllYearsPlot'}
                    % Creates a Lineplot of different variables and bar of the clibration variables
                    configs.savename = GeneralPlotting(configs.SingleRunAnalysis,'Dont Save',data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,r,1);

                case {'Regression and Scatter'}
                    % Creates the scatter plots and regression seen in Girod et al. 2009
                    configs.savename = GeneralPlotting(configs.SingleRunAnalysis,'Dont Save',data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,r,1);

                case {'Sensitivity Analysis'}
                    configs.savename = GeneralPlotting(configs.SingleRunAnalysis,'Dont Save',data{run}.cfg,exvars,data{run}.depvars,data{run}.indepvars,r,1);
                    % Post-calibration sensitivity analysis

                case {'Plot input data'}
                    i=1; fields=fieldnames(exvars); figure(1); 
                    for n=1:length(fields); for m=1:size(exvars.(fields{n}),3);
                        subplot(7,7,i); plot(squeeze(exvars.(fields{n})(r,:,m))'); title([fields{n} ' ' num2str(m)]); i=i+1;
                    end; end
            end
            savegraphs(configs)
        end
    end
end
    
if ~strcmp(configs.MultiRunAnalysis,'None')% Used to skip this analysis tools
%% Multiple run analyses, first combine the data
    
if configs.AllRegions == 1; configs.Regions = cfg.Regions; end

fields = fieldnames(data{1});
fields = fields (~strcmp(fields,'cfg')&~strcmp(fields,'exvars')&~strcmp(fields,'errors'));
for i = 1:length(fields)
    eval(['varfields = fieldnames(data{1}.' fields{i} ');']);
    varfields = varfields(~strcmp(varfields,'errors'));
    for j = 1:length(varfields)
        disp(['Merging data for ' fields{i} '.' varfields{j}]);
        eval([fields{i} '.' varfields{j} '= zeros([length(configs.Runs) size(data{1}.(fields{i}).(varfields{j})) ]);']);
        for k = 1:length(configs.Runs)
            eval([fields{i} '.' varfields{j} '(k,:,:,:) = data{configs.Runs(1)-1+k}.(fields{i}).(varfields{j});']);
        end
    end
end

varfields = fieldnames(data{1}.depvars.errors);
    for j = 1:length(varfields)
        disp(['Merging data for depvars.errors.' varfields{j}]);
        eval(['depvars.errors.' varfields{j} '= zeros([length(configs.Runs) size(data{1}.depvars.errors.(varfields{j})) ]);']);
        for k = 1:length(configs.Runs)
            eval(['depvars.errors.' varfields{j} '(k,:,:,:) = data{configs.Runs(1)-1+k}.depvars.errors.(varfields{j});']);
        end
end

% rank the data
%% Outlier detection
if configs.Removeoutliers
 [y ] = (depvars.TravelDemand(:,:,end,:));
 q = repmat(permute(quantile(y,[0.05,0.25, 0.5 0.75 0.95],1),[5 2 3 4 1]),[size(depvars.TravelDemand,1) 1 1 1 1]);
 IQR = q(:,:,:,:,4)-q(:,:,:,:,2);
 ol = y<q(:,:,:,:,2)-1.5*IQR;
 ou = y>q(:,:,:,:,4)+1.5*IQR;
 configs.outliers = ~(sum(ol,4)>0 | sum(ou,4)>0)|~(sum(y<q(:,:,:,:,1),4)>0 | sum(y>q(:,:,:,:,5),4)>0);
else
 configs.outliers = logical(ones(size(1)));
end
 
%% Next allow the user to select a datarange
if ~isfield(cfg,'MonteCarloFields'); cfg.MonteCarloFields={'CostDistrParameter', 'PriceElasticity','IncomeElasticity'}; end
r = configs.Regions(1);
configs.mcnumel   = cellfun(@(x) size(indepvars.(x),4),cfg.MonteCarloFields);
configs.variables = {};
configs.mccount   = [];
configs.tranges   = [];
for n = 1 :length(configs.mcnumel)
    configs.mccount(end+1:end+configs.mcnumel(n))   = 1:configs.mcnumel(n);
    configs.variables(end+1:end+configs.mcnumel(n)) = cfg.MonteCarloFields(n);
    configs.tranges(end+1:end+configs.mcnumel(n),1:2)= ...
        [permute(min(indepvars.(cfg.MonteCarloFields{n})(:,r,1,:),[],1)/2 ...
        +max(indepvars.(cfg.MonteCarloFields{n})(:,r,1,:),[],1)/2, [4 1 2 3]) ...
        permute((max(indepvars.(cfg.MonteCarloFields{n})(:,r,1,:),[],1)...
        -min(indepvars.(cfg.MonteCarloFields{n})(:,r,1,:),[],1))/2,[4 1 2 3])];
end

configs.variablesnum = {};
for z = 1:length(configs.variables)
configs.variablesnum{z} = [configs.variables{z},  num2str(configs.mccount(z))] ;
end

configs.ranges    = NumberEntryMenu(configs.variablesnum ,num2cell(configs.tranges,2));
configs.varlogical=logical(ones([size(indepvars.(configs.variables{1}),1),size(indepvars.(configs.variables{1}),2)]));
z=1;m=1;
for n = 1:length(configs.variables)
    if size(indepvars.(configs.variables{n}),4)>1 && m<size(indepvars.(configs.variables{n}),4);
        z = z+0.99;
    elseif size(indepvars.(configs.variables{n}),4)>1
        z=1.99;
    else
        z=1;
    end
    m=floor(z);
        configs.templogical{n} = indepvars.(configs.variables{n})(:,:,1,m) >= (configs.ranges{n}(1)-configs.ranges{n}(2))...
            & indepvars.(configs.variables{n})(:,:,1,m) <=(configs.ranges{n}(1)+configs.ranges{n}(2));
        
        configs.errorthreshold = depvars.errors.TDfit < configs.ErrorThreshold; % Filters based on errors
    
    configs.varlogical     = configs.templogical{n} & configs.errorthreshold & configs.varlogical & configs.outliers ; %& configs.varlogical
    disp([num2str(sum(configs.varlogical)) ', variables left ' configs.variables{n} ', Mode ' num2str(m)]);
end
configs.varlogical=configs.varlogical';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Batch analysis tools (requires multiple runs of data) %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if ishandle(12);set(12,'Units','centimeter'); cornerpos = get(12,'Position'); close(12); end;
        
        switch configs.MultiRunAnalysis
                case {'CO2 and Energy all regions'}
        %% CO2 and Energy intensity projections
        if isfield(data{1}.indepvars,'UseModes')
            modes               = data{1}.indepvars.UseModes(1,:);
            modes(1)            = false;
        else
            modes =[false true true true];
        end
        indepvars.EnergyIntensity         = exvars.EnergyIntensity;
        indepvars.EnergyIntensity (:,:,2) = repmat(exvars.EnergyIntensity(:,1,2),[1 130 1])...
            .*exvars.EnergyIntensity(:,:,3)./repmat(exvars.EnergyIntensity(:,1,3),[1 130 1]);
        
        indepvars.CO2Intensity         = exvars.CO2Intensity;
        indepvars.CO2Intensity (:,:,2) = repmat(exvars.CO2Intensity(:,1,2),[1 130 1])...
            .*exvars.CO2Intensity(:,:,3)./repmat(exvars.CO2Intensity(:,1,3),[1 130 1]);
        
        years = 35:130;
        
            for r=configs.Regions
                if  sum(configs.varlogical (r,:))==0
                    disp('skipping region');
                else
                    configs.savename    = [configs.MultiRunAnalysis ' CO2 ' cfg.shortregions{r} ' ' configs.OutputFile];
                    PlotRange           = 35:130;
                    fanput              = depvars.TravelDemand(:,r,years,modes).*repmat(permute(...
                        indepvars.CO2Intensity(r,years,modes),[4 1 2 3])./10^6,[size( configs.varlogical(:,:),2) 1 1 1]);
                    fanput(~configs.varlogical) = nan;
                    fanput(fanput>1000) = nan;
                    FanChart(fanput,[cfg.modes],4,cfg,r,true,true,[1/5 1/5]);
                    output.(cfg.shortregions{r}) =fanput;
                end
                
                title(['Co2 eq emissions ' cfg.shortregions{r} '(n=' num2str(sum(~isnan(fanput(:,1,1,1)))) ')']);
                ylabel('CO2 eq Emissions (t/Cap/yr)');
                savegraphs(configs)
            end
            
            [ a b] = mkdir('Output/Energyregion');
            save(['Output/Energyregion/' cfg.Model '.mat'],'output');
                
            for r=configs.Regions
                if  sum(configs.varlogical (r,:))==0
                    disp('skipping region');
                else
                configs.savename    = [configs.MultiRunAnalysis ' GJ ' cfg.shortregions{r} ' ' configs.OutputFile];
                PlotRange           = 35:130;
                fanput              = depvars.TravelDemand( :,r,years,modes).*repmat(permute...
                    (indepvars.EnergyIntensity(r,years,modes),[4 1 2 3])./10^3,[size( configs.varlogical(:,:),2) 1 1 1]);
                fanput(fanput>1000) = nan;
                fanput(~configs.varlogical) = nan;
                FanChart(fanput,[cfg.modes],4,cfg,r,true,true,[1/5 1/5]);
                output.(cfg.shortregions{r}) =fanput;
                
                title(['Energy use ' cfg.shortregions{r} '(n=' num2str(sum(~isnan(fanput(:,1,1,1)))) ')']);
                ylabel('Energy use (GJ/Cap/yr)');
                savegraphs(configs)
                end
            end
            
            [ a b] = mkdir('Output/CO2region');
            save(['Output/CO2region/' cfg.Model '.mat'],'output');
            
            for r=configs.Regions
                if  sum(configs.varlogical (r,:))==0
                    disp('skipping region');
                else
                configs.savename    = [configs.MultiRunAnalysis ' TD ' cfg.shortregions{r} ' ' configs.OutputFile];
                PlotRange           = 35:130;
                fanput              = depvars.TravelDemand( :,r,years,modes);
                fanput(fanput>100000) = nan;                    
                fanput(~configs.varlogical) = nan;

                FanChart(fanput,[cfg.modes],4,cfg,r,true,true,[1/5 1/5]);
                output.(cfg.shortregions{r}) =fanput;
                
                title(['Travel demand ' cfg.shortregions{r} '(n=' num2str(sum(~isnan(fanput(:,1,1,1)))) ')']);
                ylabel('Travle demand (pkm/yr)');
                savegraphs(configs)
                end
            end
            [ a b] = mkdir('Output/TDregion');
            save(['Output/TDregion/' cfg.Model '.mat'],'output');

        case {'CO2 and Energy intensity projections'}
 %%%% Repeat for all regions combines %%%%%%%
 
        if isfield(data{1}.indepvars,'UseModes')
            modes               = data{1}.indepvars.UseModes(1,:);
            modes(1)            = false;
        else
            modes =[false true true true];
        end
        
        indepvars.EnergyIntensity         = exvars.EnergyIntensity;
        indepvars.EnergyIntensity (:,:,2) = repmat(exvars.EnergyIntensity(:,1,2),[1 130 1])...
            .*exvars.EnergyIntensity(:,:,3)./repmat(exvars.EnergyIntensity(:,1,3),[1 130 1]);
        
        indepvars.CO2Intensity         = exvars.CO2Intensity;
        indepvars.CO2Intensity (:,:,2) = repmat(exvars.CO2Intensity(:,1,2),[1 130 1])...
            .*exvars.CO2Intensity(:,:,3)./repmat(exvars.CO2Intensity(:,1,3),[1 130 1]);
        
        years = 35:130;
        
        for n=1:4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CO2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               cfg.shortregions  = {'Industrialized','Reforming','Developing','World'};
               temp            = ~logical(1:12);
               temp2           = temp;
               cfg.regionrange = {[1 2 3],[4 5],[6:11],1:11};
               temp(cfg.regionrange{n}) = true;
               temp2(configs.Regions)   = true;
               r  = find(temp & temp2);
               
                configs.savename    = [configs.MultiRunAnalysis ' CO2 ' cfg.shortregions{n} '-' configs.OutputFile];
                PlotRange           = 35:130;
                fanput              = sum((depvars.TravelDemand(:,r,years,modes).*...
                    repmat(permute(indepvars.CO2Intensity(r,years,modes),[4 1 2 3])./10^6,[size(depvars.TravelDemand,1) 1 1 1])).*...
                    repmat(permute(exvars.Population(r,years)*10^6,[3 1 2]),[size(depvars.TravelDemand,1) 1 1 sum(modes)]),2)...
                    ./sum(repmat(permute(exvars.Population(r,years)*10^6,[3 1 2]),[size(depvars.TravelDemand,1) 1 1 sum(modes)]),2);
                fanput(repmat(sum(configs.varlogical(r,:)',2)~=length(r),[1 1 size(fanput,3) size(fanput,4)])) = nan;
                fanput(fanput>1000) = nan;
                FanChart(fanput,[cfg.modes],4,cfg,n,true,true,[1/5 1/5]);
                
                title(['Co2 eq em ' cfg.shortregions{n} '(n=' num2str(sum(~isnan(fanput(:,1,1,1)))) ')']);
                ylabel('CO2 eq em (t/Cap/yr)');
                savegraphs(configs)
                                output.(cfg.shortregions{n}) = fanput;

        end
        
            save(['Output/CO2/' cfg.Model '.mat'],'output');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ENERGY%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for n = 1:4
                
               cfg.shortregions  = {'Industrialized','Reforming','Developing','World'};
               temp            = ~logical(1:12);
               temp2           = temp;
               cfg.regionrange = {[1 2 3],[4 5],[6:11],1:11};
               temp(cfg.regionrange{n}) = true;
               temp2(configs.Regions)   = true;
               r  = find(temp & temp2);
                configs.savename    = [configs.MultiRunAnalysis ' GJ ' cfg.shortregions{n} ' ' configs.OutputFile];
                PlotRange           = 35:130;
                fanput              = sum((depvars.TravelDemand(:,r,years,modes).*...
                    repmat(permute(indepvars.EnergyIntensity(r,years,modes),[4 1 2 3])./10^6,[size(depvars.TravelDemand,1) 1 1 1])).*...
                    repmat(permute(exvars.Population(r,years)*10^6,[3 1 2]),[size(depvars.TravelDemand,1) 1 1 sum(modes)]),2)...
                    ./sum(repmat(permute(exvars.Population(r,years)*10^6,[3 1 2]),[size(depvars.TravelDemand,1) 1 1 sum(modes)]),2);
                fanput(repmat(sum(configs.varlogical(r,:)',2)~=length(r),[1 1 size(fanput,3) size(fanput,4)])) = nan;
                fanput(fanput>1000) = nan;
                 FanChart(fanput,[cfg.modes],4,cfg,n,true,true,[1/5 1/5]);
                               
                title(['Energy use ' cfg.shortregions{n} '(n=' num2str(sum(~isnan(fanput(:,1,1,1)))) ')']);
                ylabel('Energy use (GJ/Cap/yr)');
                savegraphs(configs)
                output.(cfg.shortregions{n}) = fanput;
        end
                    save(['Output/Energy/' cfg.Model '.mat'],'output');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%TRAVELDEMAND%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
                    for n = 1:4
                
               cfg.shortregions  = {'Industrialized','Reforming','Developing','World'};
               temp            = ~logical(1:12);
               temp2           = temp;
               cfg.regionrange = {[1 2 3],[4 5],[6:11],1:11};
               temp(cfg.regionrange{n}) = true;
               temp2(configs.Regions)   = true;
               r  = find(temp & temp2);
                configs.savename    = [configs.MultiRunAnalysis ' TD ' cfg.shortregions{n} ' ' configs.OutputFile];
                PlotRange           = 35:130;
                fanput              = sum(depvars.TravelDemand(:,r,years,modes)./1000.*repmat(permute(exvars.Population(r,years),[3 1 2]),[size(depvars.TravelDemand,1) 1 1 sum(modes)]),2)...
                    ./sum(repmat(permute(exvars.Population(r,years),[3 1 2]),[size(depvars.TravelDemand,1) 1 1 sum(modes)]),2);
                fanput(repmat(sum(configs.varlogical(r,:)',2)~=length(r),[1 1 size(fanput,3) size(fanput,4)])) = nan;
                fanput(fanput>100) = nan;
                 FanChart(fanput,[cfg.modes],4,cfg,n,true,true,[1/5 1/5]);
                               
                title(['Travel demand ' cfg.shortregions{n} '(n=' num2str(sum(~isnan(fanput(:,1,1,1)))) ')']);
                ylabel('TravelDemand (thousand pkm/yr)');
                savegraphs(configs)
                output.(cfg.shortregions{n}) = fanput;
                    end
            [a b] = mkdir('Output/TD');
            save(['Output/TD/' cfg.Model '.mat'],'output');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%END%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           
            case {'Boxplot indepvars per region and mode'}
                aboxplot(permute(a,[4 1 2 3]),'Colorgrad','green_up')
                a=1/(indepvars.PreferenceFactorStart(:,:,:,:).^0.3./(repmat(indepvars.PreferenceFactorStart(:,:,:,3),[1 1 1 4])).^0.3);
                for n = [ 1 2 4]; subplot(3,1,ceil(n/1.5)); aboxplot(permute(1./a(:,:,:,n),[4 1 2 3]),'Colorgrad','green_up');
                    title(['Preference Factor for ' cfg.modes{n}]); set(gca,'Ylim',[m,5],'YTick',0:5,'Xticklabel',cfg.shortregions); ylabel('Preference Factor'); end
                
            case {'BoxPlots of NRSME'}
                errorfields = fieldnames(depvars.errors);
                selection = SelectionMenu(fieldnames(depvars.errors),'Select error indicators to compare');
                selectioncell = errorfields(selection);
                for n=1:length(selectioncell)
                    globals.fid = figure(12);
                    cfg.shortregions = {'NAM','PAO','WEU','EEU','FSU','Afr','CPA','LAM','MEA','PAS','SA','Wor'};
                    aboxplot(depvars.errors.(selectioncell{n}),'WidthS',1,'Colorgrad','green_up','Colorrev',true,'Labels',cfg.shortregions);
                    ylabel(['NRSME of ' selectioncell{n}]); xlabel('Region'); title(['Boxplot of error distribution for Monte Carlo simulation; ' selectioncell{n}]);
                    pubgraph(globals.fid,12,1.5,'white'); set(gca,'YLim',[0 0.6]);
                    configs.savename = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' )];
                    savegraphs(configs)
                end
                
            case {'BoxPlot of indepvar per mode'}
            %% Boxplot over modes
            for r=configs.Regions
                
            xfields                 = fieldnames(indepvars);
            selection               = SelectionMenu(fieldnames(indepvars),'Select error indicators to compare');
            cfg.XFields             = xfields(selection);

                for x = 1:length(cfg.XFields)
                        globals.fid = figure(12);
                        corrcoef(permute(indepvars.(cfg.XFields{1})(:,12,1,:),[1 4 2 3]))
                        aboxplot(permute(indepvars.(cfg.XFields{x})(:,r,1,:),[1 4 2 3]),'WidthS',1,'Colorgrad','green_up','Colorrev',true,'Labels',cfg.modes);
                        ylabel([cfg.XFields{x}]); 
                        xlabel('Region'); title(['Distribution of ' cfg.XFields{x} ', ' cfg.allregions{r}]);
                        pubgraph(globals.fid,12,1.5,'white');
                         set(gca,'YLim',[min(min(min(min(indepvars.(cfg.XFields{1})(:,:,:,:))))) max(max(max(max(indepvars.(cfg.XFields{1})(:,:,:,:)))))]);
                        configs.savename = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' boxplot of ' cfg.XFields{x}];
                        savegraphs(configs)
                        close (12);
                end
            end
                            
            case {'Send to workspace'}
                sendtoworkspace(cfg,exvars,depvars,indepvars);
                
            case {'BoxPlot of indepvars'}
            xfields                 = fieldnames(indepvars);
            selection               = SelectionMenu(fieldnames(indepvars),'Select error indicators to compare');
            cfg.XFields             = xfields(selection);
            
            for x = 1:length(cfg.XFields)
                for y = 1:size(indepvars.(cfg.XFields{x}),4)
                    globals.fid = figure(12);
                    subplot(size(indepvars.(cfg.XFields{x}),4),1,y);
                    cfg.shortregions = {'NAM','PAO','WEU','EEU','FSU','Afr','CPA','LAM','MEA','PAS','SA','Wor'};
                    aboxplot(indepvars.(cfg.XFields{x})(configs.varlogical(r,:),:,1,y),'WidthS',1,'Colorgrad','green_up','Colorrev',true,'Labels',cfg.shortregions);
                    ylabel([cfg.XFields{x} num2str(y)]); xlabel('Region'); title(['Boxplot of distribution for Monte Carlo simulation; ' cfg.XFields{x} num2str(y)]);
                    pubgraph(globals.fid,12,1.5,'white');
                    set(gca,'YLim',[min(min(min(min(indepvars.(cfg.XFields{1})(:,:,:,:))))) max(max(max(max(indepvars.(cfg.XFields{1})(:,:,:,:)))))]);
                    configs.savename = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' NRSME of ' cfg.XFields{x} num2str(y)];
                   savegraphs(configs)
                end
            end
                
            case {'Scatter Plots, errors added'}
               
            %% Generates some simple scatter plots
            errorfields = fieldnames(depvars.errors);
            selection = SelectionMenu(fieldnames(depvars.errors),'Select error indicators to compare');
            selectioncell = errorfields(selection);
            depvars.errors.error = ones([size(depvars.errors.(selectioncell{1}))]);
            configs.errorsadded   = '';
            for n=1:length(selectioncell)
                depvars.errors.error = depvars.errors.(selectioncell{n}) .* depvars.errors.error;
                configs.errorsadded = [configs.errorsadded selectioncell{n} ' * '];
            end
            configs.errorsadded = configs.errorsadded(1:end-3);
                for r = configs.Regions

                    t           = 1 ;
                    cfg.MonteCarloFields    = {'CostDistrParameter', 'PriceElasticity','IncomeElasticity'}; % REMOVE
                    for x= 1:length(cfg.MonteCarloFields)
                        globals.fid = figure(12);
                        scatter(indepvars.(cfg.MonteCarloFields{x})(:,r),depvars.errors.error(:,r),'x')
                        xlabel (cfg.MonteCarloFields{x});
                        ylabel(['NRSME of ' selectioncell{n}]); 
                        title(['Region ' cfg.allregions{r} ', NRSME of ' configs.errorsadded])
                        pubgraph(globals.fid,14,1.5,'white')
                        axis([min(indepvars.(cfg.MonteCarloFields{x})(:,r)) max(indepvars.(cfg.MonteCarloFields{x})(:,r)) 0 max(depvars.errors.error(:,r))])
                        configs.savename = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ', Variable ' cfg.MonteCarloFields{x} ', Region ' cfg.allregions{r}];
                        savegraphs(configs)
                    end
                 end
                
                
            case {'Scatter Plots, independant errors'}
               
            %% Generates some simple scatter plots
            errorfields = fieldnames(depvars.errors);
            selection = SelectionMenu(fieldnames(depvars.errors),'Select error indicators to compare');
            selectioncell = errorfields(selection);
            for n=1:length(selectioncell)
                for r = configs.Regions

                    t           = 1 ;
      
                    for x= 1:length(unique(cfg.MonteCarloFields))
                        for yy=1:size(indepvars.(cfg.MonteCarloFields{x}),4);
                        globals.fid = figure(12);
                        scatter(indepvars.(cfg.MonteCarloFields{x})(configs.varlogical(r,:),r,1,yy),depvars.errors.(selectioncell{n})(configs.varlogical(r,:),r),'x','MarkerFaceColor',[0.2 0.2 0.2],'MarkerEdgeColor',[0.2 0.2 0.2])
                        box on
                        xlabel (cfg.MonteCarloFields{x});
                        ylabel(['NRSME of ' selectioncell{n}]); 
                        if size(indepvars.(cfg.MonteCarloFields{x}),4) == 1
                        title(['Region ' cfg.allregions{r} ', NRSME of ' selectioncell{n}])
                        else
                        title(['Region ' cfg.allregions{r} ' Mode , ' cfg.modes{yy} ', NRSME of ' selectioncell{n}])
                        end
                        pubgraph(globals.fid,14,1.5,'white')
                        axis([min(indepvars.(cfg.MonteCarloFields{x})(:,r,1,yy)) max(indepvars.(cfg.MonteCarloFields{x})(:,r,1,yy))+0.001 -0.001 max(depvars.errors.(selectioncell{n})(:,r))])
                        configs.savename = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ', Variable ' cfg.MonteCarloFields{x} ', Region ' cfg.allregions{r} num2str(yy)];
                        savegraphs(configs)
                        end
                    end
                 end
            end
            
            case {'Scatter plots, select variables'}
            %% Generates some simple scatter plots
            errorfields             = fieldnames(depvars.errors);
            selection               = SelectionMenu(fieldnames(depvars.errors),'Select error indicators to compare');
            selectioncell           = errorfields(selection);
            xfields                 = fieldnames(indepvars);
            selection               = SelectionMenu(fieldnames(indepvars),'Select error indicators to compare');
            cfg.XFields             = xfields(selection);
            for n=1:length(selectioncell)
                for r = configs.Regions
                    t           = 1 ;
                    for x= 1:length(cfg.XFields)
                        for z = 1:size(indepvars.(cfg.XFields{x}),4)
                        globals.fid = figure(12);
                        scatter(indepvars.(cfg.XFields{x})(:,r,1,z),depvars.errors.(selectioncell{n})(:,r),'x')
                        xlabel ([cfg.XFields{x} num2str(z)]);
                        ylabel(['NRSME of ' selectioncell{n}]); 
                        title(['Region ' cfg.allregions{r} ', NRSME of ' selectioncell{n}])
                        pubgraph(globals.fid,14,1.5,'white')
                        axis([min(indepvars.(cfg.XFields{x})(:,r,1,z)) max(indepvars.(cfg.XFields{x})(:,r,1,z)) 0 max(depvars.errors.(selectioncell{n})(:,r))])
                        configs.savename = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ', Variable ' cfg.XFields{x} num2str(z) ', Region ' cfg.allregions{r}];
                        savegraphs(configs)
                        end
                    end
                end
            end
                
            case {'Probablistic Projection'}
                %% Probablistic projections, first take only valid data                 
                    % Does the actual plotting
                for r = configs.Regions
                    cfg.ShowLegend = true;
                    if sum(configs.varlogical(r,:)) >1
                    configs.savename = FanChart(depvars.TravelDemand(configs.varlogical(r,:),r,35:130,:)/1000,cfg.modes,4,cfg,r);
                    savegraphs(configs)
                    else
                        disp(['Region ' cfg.allregions{r} ' doesnt have any valid points']);
                    end
                    
                end
        end
end
        %% Batch over multiple files
else
    load('C:\Users\Roel\Documents\Master thesis ETH\Data\exvars 1971-2100 Original Dataset Regions11Modes4.mat')
    
    options      = {'TD','CO2','Energy'};
    investigated   = options(SelectionMenu(options));
    options        = {'yes','no'};
    multipopulation=  options(SelectionMenu(options,'Multiply by population?'));
    for invs = 1:length(investigated)
    switch multipopulation{1}
        case 'yes'
            
            population    = {sum(exvars.Population(1:3,35:130),1)*10^6,sum(exvars.Population(4:5,35:130),1) * 10^6 ...
                ,sum(exvars.Population(6:11,35:130),1) *10^6,sum(exvars.Population(1:11,35:130),1)*10^6};
            switch investigated{invs}
                case 'TD'
                    plottitle = ['Total Travel demand ' ];
                    labely    = ['Travel demand (Tera pkm/yr)'];
                    multiplier= 1/10^12;
                case 'CO2'
                    plottitle = ['Total CO2 eq emissions ' ];
                    labely    = ['Total CO2 eq emissions (Gt/yr)'];
                    multiplier= 1/10^9;
                case 'Energy'            
                    plottitle = ['Total Energy use ' ];
                    multiplier= 1/10^9;
                    labely    = ['Energy Use (EJ/yr)'];
            end
            
        case 'no'
                        population     = {sum(exvars.Population(1:3,35:130),1),sum(exvars.Population(4:5,35:130),1)...
                ,sum(exvars.Population(6:11,35:130),1),sum(exvars.Population(1:11,35:130),1)};
            switch investigated{invs}
                case 'TD'
                    plottitle = ['Total Travel demand' ];
                    labely    = ['Travel demand (pkm/cap/yr)'];
                    multiplier= 1000;
                case 'CO2'
                    plottitle = ['Total CO2 emissions' ];
                    labely    = ['Total CO2 eq emissions (t/cap/yr)'];
                    multiplier= 1;
                case 'Energy'            
                    plottitle = ['Total Energy use' ];
                    multiplier= 1;
                    labely    = ['Energy Use (GJ/cap/yr)'];
            end
    end
    
    switch configs.MultiFileAnalysis
        case {'Boxplot all regions'}
            if strcmp(multipopulation,'yes')
            population = num2cell(exvars.Population([1:6 8:11],35:130).*10^6,2); else
            population = num2cell(exvars.Population([1:6 8:11],35:130),2); end
        cfg.files      = {'ModelTravelLinearTW' 'ModelGcam' 'ModelSimple' 'ModelPoles' };
            for file = 1:4
                output = load(['Output/' investigated{invs} 'region/' cfg.files{file} '.mat']);
                Merged{file} = output.output;
                disp (['loaded ' cfg.files{file}]);
                clear output;
            end
             allregions = {'North America','Pacific OECD','Western Europe','Eastern Europe',...
                    'Former Soviet Union','Sub-Saharan Africa'... %,'Centrally Planned Asia',...
                    'Latin America','Middle East & North Africa','Other Pacific Asia',...
                    'South Asia','World'};
            cfg.modes      = {'Timer-Travel' 'Gcam' 'Simple' 'Poles'};
            fields = fieldnames(Merged{1});
            toplot = logical([[1 1 1 1];[1 1 1 1];[1 1 1 1];...
                [1 1 1 0]; [1 1 1 0]; ...
                [1 1 0 0];[1 1 0 0];[1 1 0 0];[1 1 0 0];[1 1 0 0];]);
            for n = 1:length(fields)
                fanput.(fields{n}) = zeros([5000 1 size(Merged{1}.(fields{n}),3) length(Merged)]);
                for m=1:length(Merged)
                    fanput.(fields{n})(:,:,:,m) = repmat(sum(Merged{m}.(fields{n}),4),[5000/size(Merged{m}.(fields{n}),1) 1 1 1])...
                        .* permute(repmat(population{n},[5000 1 1]),[1 3 2])*multiplier;
                end
                cfg.allregions = fields;
                cfg.Time       = 1:130;
                cfg.years      = 1971:2100;
                PlotRange      = 1:130;
                cfg.ShowLegend=false;
                FanChart(fanput.(fields{n})(:,:,:,toplot(n,:)),[cfg.modes(toplot(n,:))],4,cfg,n,true,true,[2/5 0.75/5]);
                ylabel(labely);

                title([allregions{n}]);
                 configs.savename = [investigated{invs} ' production summarized ' fields{n} datestr(now,'yyyy-mm-dd HH-MM-SS')];
                     savegraphs(configs)
            end  
            
        case {'Probablistic 4 regions'}   
                cfg.files      = {'ModelTravelLinearTW' 'ModelGcam' 'ModelSimple' 'ModelPoles' };
    for file = 1:4
        output = load(['Output/' investigated{invs} '/' cfg.files{file} '.mat']);
        Merged{file} = output.output;
        disp (['loaded ' cfg.files{file}]);
        clear output;
    end
            cfg.modes      = {'Timer-Travel' 'Gcam' 'Simple' 'Poles'};
            fields = fieldnames(Merged{1});
            toplot = logical([[1 1 1 1]; [1 1 1 0]; [1 1 0 0]; [1 1 0 0]]);
            for n = 1:length(fields)
                fanput.(fields{n}) = zeros([5000 1 size(Merged{1}.(fields{n}),3) length(Merged)]);
                for m=1:length(Merged)
                    fanput.(fields{n})(:,:,:,m) = repmat(sum(Merged{m}.(fields{n}),4),[5000/size(Merged{m}.(fields{n}),1) 1 1 1])...
                        .* permute(repmat(population{n},[5000 1 1]),[1 3 2])*multiplier;
                end
                cfg.allregions = fields;
                cfg.Time       = 1:130;
                cfg.years      = 1971:2100;
                cfg.ShowLegend=true;
                FanChart(fanput.(fields{n})(:,:,:,toplot(n,:)),[cfg.modes(toplot(n,:))],4,cfg,n,true,true,[2/5 0.75/5]);
                ylabel(labely);
                title([plottitle fields{n}]);
                 configs.savename = [investigated{invs} ' production summarized ' fields{n} datestr(now,'yyyy-mm-dd HH-MM-SS')];
                     savegraphs(configs)
            end
            end
        end
    end
end

function sendtoworkspace(cfg,exvars,depvars,indepvars)
%% Exports the variables to workspace, for one off graphical magic

assignin('base', 'cfg', cfg);
assignin('base', 'exvars', exvars);
assignin('base', 'depvars', depvars);
assignin('base', 'indepvars', indepvars);

end

function savegraphs(configs)
%% Save graphs
global globals
set(globals.fid,'Units','centimeter');
configs.savename = [configs.savename configs.OutputFile];
if strcmp('Default',configs.PlotExportType)==1
        disp(['Saving graph as ' configs.savename]);
        [trash ~] = mkdir(['Output/' configs.OutputFolder]);
        saveas(globals.fid,['Output/' configs.OutputFolder '/' configs.savename '.jpg']);

else
    
    %% Suitable for publication
    configs.fulldir = pwd;
    switch configs.PlotExportType    
        case 'A4'
                pubgraph(globals.fid,10,1.5,'white') % makes the graphs a little more readable
                set(globals.fid,'Units','centimeter')
                set(globals.fid,'Position',get(globals.fid,'PaperPosition').*[0 0 1 1]+get(globals.fid,'Position').*[1 1 0 0]);
                
        case 'none'
                
        case 'Publication Tiny'
                pubgraph(globals.fid,10,1,'white') % makes the graphs a little more readable
                set(globals.fid,'Units','centimeter')
                set(globals.fid,'Position',[0 0 5.5 5.5]+get(globals.fid,'Position').*[1 1 0 0]);
                set(globals.fid,'PaperPosition',[0 0 5.5 5.5]);
                
        case 'Publication Small Flat'
                pubgraph(globals.fid,10,1,'white') % makes the graphs a little more readable
                set(globals.fid,'Units','centimeter')
                set(globals.fid,'Position',[0 0 8 6]+get(globals.fid,'Position').*[1 1 0 0]);
                set(globals.fid,'PaperPosition',[0 0 8 6]);
                
        case 'Publication Small'
                pubgraph(globals.fid,10,1.5,'white') % makes the graphs a little more readable
                set(globals.fid,'Units','centimeter')
                set(globals.fid,'Position',[0 0 8 8]+get(globals.fid,'Position').*[1 1 0 0]);
                set(globals.fid,'PaperPosition',[0 0 8 8]);

        case 'Publication Large' 
                pubgraph(globals.fid,10,1.5,'white') % makes the graphs a little more readable
                set(globals.fid,'Units','centimeter')
                set(globals.fid,'Position',[0 0 16 8]+get(globals.fid,'Position').*[1 1 0 0]);
                set(globals.fid,'PaperPosition',[0 0 16 8]);
    end
    
    [a b] = mkdir([configs.fulldir '/Publishable/' configs.OutputFolder]);
    print('-depsc','-tiff','-r300',['-f' num2str(globals.fid)],[configs.fulldir '/Publishable/' configs.OutputFolder '/' configs.savename '.eps']);
% cd('Export figure');
%     export_fig([configs.fulldir '/Publishable/' configs.OutputFolder '/' configs.savename '-ef.eps'],'-transparent');
% cd ..
    disp(['Saved as ' configs.savename '.eps']);
end

end

function aboxplot(X,varargin)

% Parameters
widthl = 0.7;
widths = 0.8;
widthe = 0.4;
outMarker = '.';
outMarkerSize = 3;
outMarkerEdgeColor = [0.6 0.6 0.6];
outMarkerFaceColor = [0.6 0.6 0.6];
alpha = 0.05;
cmap = [];
colorrev = 0;
colorgrd = 'blue_down';

% Get the number or data matrices
if iscell(X)
    d = length(X);
else
    % If data is a matrix extend to a 3D array
    if 2 == ndims(X)
        X = reshape(X, [1,size(X)]);
    end
    d = size(X,1);
end;

% Get the data size
if iscell(X)
    n = size(X{1},2);
else
    n = size(X,3);
end

% Set the labels
labels = cell(n,1);
for i=1:n
    labels{i} = num2str(i);
end

% Optional arguments
optargin = size(varargin,2);

i = 1;
while i <= optargin
    switch lower(varargin{i})
        case 'labels'
            labels = varargin{i+1};
        case 'colormap'
            cmap = varargin{i+1};
        case 'colorgrad'
            colorgrd = varargin{i+1};
        case 'colorrev'
            colorrev = varargin{i+1};
        case 'outliermarker'
            outMarker = varargin{i+1};
        case 'outliermarkersize'
            outMarkerSize = varargin{i+1};
        case 'outliermarkeredgecolor'
            outMarkerEdgeColor = varargin{i+1};
        case 'outliermarkerfacecolor'
            outMarkerFaceColor = varargin{i+1};
        case 'widthl'
            widthl = varargin{i+1};
        case 'widths'
            widths = varargin{i+1};
        case 'widthe'
            widthe = varargin{i+1};
    end
    i = i + 2;
end

% Colors
colors = cell(d,n);

if colorrev
    %  Set colormap
    if isempty(cmap)
        cmap = colorgrad(n,colorgrd);
    end
    if size(cmap,1) ~= n
        error('The number of colors in the colormap must equal n.');
    end
    for j=1:d
        for i=1:n
            colors{j,i} = cmap(i,:);
        end
    end
else
    %  Set colormap
    if isempty(cmap)
        cmap = colorgrad(d,colorgrd);
    end
    if size(cmap,1) ~= d
        error('The number of colors in the colormap must equal n.');
    end
    for j=1:d
        for i=1:n
            colors{j,i} = cmap(j,:);
        end
    end
end

xlim([0.5 n+0.5]);

hgg = zeros(d,1);

for j=1:d
    % Get the j matrix
    if iscell(X)
        Y = X{j};
    else
        Y = squeeze(X(j,:,:));
    end
    
    % Create a hggroup for each data set
    hgg(j) = hggroup();
    set(get(get(hgg(j),'Annotation'),'LegendInformation'),'IconDisplayStyle','on');
    legendinfo(hgg(j),'patch',...
        'LineWidth',0.5,...
        'EdgeColor','k',...
        'FaceColor',colors{j,1},...
        'LineStyle','-',...
        'XData',[0 0 1 1 0],...
        'YData',[0 1 1 0 0]);
    
    for i=1:n
    
        % Calculate the mean and confidence intervals
        [q1 q2 q3 fu fl ou ol] = quartile(Y(:,i));
        u = nanmean(Y(:,i));

        % large interval  [i - widthl/2 i + widthl/2] delta = widthl
        % medium interval start: i - widthl/2 + (j-1) * widthl / d
        % medium interval end: i - widthl/2 + j * widthl / d
        % medium interval width: widthl / d
        % medium interval middle: i-widthl/2+(2*j-1)*widthl/(2*d)
        % small interval width: widths*widthl/d
        % small interval start: i-widthl/2+(2*j-1-widths)*widthl/(2*d)
  
        % Plot outliers
        hold on;
        plot((i-widthl/2+(2*j-1)*widthl/(2*d)).*ones(size(ou)),ou,...
            'LineStyle','none',...
            'Marker',outMarker,...
            'MarkerSize',outMarkerSize,...
            'MarkerEdgeColor',outMarkerEdgeColor,...
            'MarkerFaceColor',outMarkerFaceColor,...
            'HitTest','off',...
            'Parent',hgg(j));
        plot((i-widthl/2+(2*j-1)*widthl/(2*d)).*ones(size(ol)),ol,...
            'LineStyle','none',...
            'Marker',outMarker,...
            'MarkerSize',outMarkerSize,...
            'MarkerEdgeColor',outMarkerEdgeColor,...
            'MarkerFaceColor',outMarkerFaceColor,...
            'HitTest','off',...
            'Parent',hgg(j));
        hold off;
        
        % Plot fence
        line([i-widthl/2+(2*j-1)*widthl/(2*d) i-widthl/2+(2*j-1)*widthl/(2*d)],[fu fl],...
            'Color','k','LineStyle',':','HitTest','off','Parent',hgg(j));
        line([i-widthl/2+(2*j-1-widthe)*widthl/(2*d) i-widthl/2+(2*j-1+widthe)*widthl/(2*d)],[fu fu],...
            'Color','k','HitTest','off','Parent',hgg(j));
        line([i-widthl/2+(2*j-1-widthe)*widthl/(2*d) i-widthl/2+(2*j-1+widthe)*widthl/(2*d)],[fl fl],...
            'Color','k','HitTest','off','Parent',hgg(j));
        
        % Plot quantile
        if q3 > q1
            rectangle('Position',[i-widthl/2+(2*j-1-widths)*widthl/(2*d) q1 widths*widthl/d q3-q1],...
                'EdgeColor','k','FaceColor',colors{j,i},'HitTest','off','Parent',hgg(j));
        end
        
        % Plot median
        line([i-widthl/2+(2*j-1-widths)*widthl/(2*d) i-widthl/2+(2*j-1+widths)*widthl/(2*d)],[q2 q2],...
            'Color','k','LineWidth',1,'HitTest','off','Parent',hgg(j));
        
        % Plot mean
        hold on;
        plot(i-widthl/2+(2*j-1)*widthl/(2*d), u,...
            'LineStyle','none',...
            'Marker','o',...
            'MarkerEdgeColor','k',...
            'MarkerFaceColor',colors{j,i},...
            'HitTest','off','Parent',hgg(j));
        hold off;
    end
end

box on;

set(gca,'XTick',1:n);
set(gca,'XTickLabel',labels);

end

function c = colorgrad(varargin)

n = 16;
t = 'blue_down';

switch length(varargin)
    case 1
        n = varargin{1};
    case 2
        n = varargin{1};
        t = varargin{2};
end

switch lower(t)
    case 'blue_up'
        c = cat(2,linspace(0,0.6,n)',linspace(0.2,0.8,n)',linspace(0.6,1,n)');
    case 'blue_down'
        c = cat(2,linspace(0.6,0,n)',linspace(0.8,0.2,n)',linspace(1,0.6,n)');
    case 'orange_up'
        c = cat(2,linspace(1,248/255,n)',linspace(0.6,224/255,n)',linspace(0,124/255,n)');
    case 'orange_down'
        c = cat(2,linspace(248/255,1,n)',linspace(224/255,0.6,n)',linspace(124/255,0,n)');
    case 'green_up'
        c = cat(2,linspace(0.2,0.6,n)',linspace(0.6,1,n)',linspace(0.2,0.6,n)');
    case 'green_down'
        c = cat(2,linspace(0.6,0.2,n)',linspace(1,0.6,n)',linspace(0.6,0.2,n)');
    case 'red_up'
        c = cat(2,linspace(.8,1,n)',linspace(.2,.6,n)',linspace(.2,.6,n)');
    case 'red_down'
        c = cat(2,linspace(1,.8,n)',linspace(.6,.2,n)',linspace(.6,.2,n)');
    otherwise
        error('No such color gradient.');
end

end


function [q1 q2 q3 fu fl ou ol] = quartile(x)

% rank the data
y = sort(x);

% compute 50th percentile (second quartile)
q2 = nanmedian(y);

% compute 25th percentile (first quartile)
q1 = nanmedian(y(y<=q2));

% compute 75th percentile (third quartile)
q3 = nanmedian(y(y>=q2));

% compute Interquartile Range (IQR)
IQR = q3-q1;

fl = min(y(y>=q1-1.5*IQR));
fu = max(y(y<=q3+1.5*IQR));

ol = y(y<q1-1.5*IQR);
ou = y(y>q3+1.5*IQR);

end